/*
 * ToroidalPore_test.cpp
 *
 *  Created on: Nov 15, 2011
 *      Author: yi
 */

#include "mbuild.h"
#include "time.h"
#include <iostream>
int testSetupPore(ToroidalPore& pore) {
	double R0,K0,T;
	double xmin,xmax,ymin,ymax,zmin,zmax;
	R0=0;K0=0;T=27;
//	cout<<"Toroidal pore size(set R0=0 and K0=0 if don't want to have a toroidal pore).\n";
//	cout<<"R0:";
//	cin>>R0;
//	cout<<"K0:";
//	cin>>K0;
//	cout<<"Membrane Thickness(T):";
//	cin>>T;
	try {
		pore.setup(R0,K0,T);
	} catch(ParameterError& error) {
      cout<<"Incorrect Parameters!\n";
      return 0;
	}
	xmin=ymin=zmin=-40;xmax=ymax=zmax=40;

//	cout<<"Xmin:";
//	cin>>xmin;
//	cout<<"Xmax:";
//	cin>>xmax;
//	cout<<"Ymin:";
//	cin>>ymin;
//	cout<<"Ymax:";
//	cin>>ymax;
//	cout<<"Zmin:";
//	cin>>zmin;
//	cout<<"Zmax:";
//	cin>>zmax;
	//Do not use fast match.
    pore.set_region(xmin,xmax,ymin,ymax,zmin,zmax);
    double cutoff=10;
    int ncyc=3;
    //TODO : setup calculation.

//    cout<<"Integration cutoff(for improving the accuracy of calculating density around the surface):";
//    cin>>cutoff;
//    cout<<"Integration cycles(more cycles,higher accuracy but longer calculation time):";
//    cin>>ncyc;
    try {
    pore.set_cutoff(cutoff);
    pore.set_integration_cyc(ncyc);
    } catch(ParameterError& error) {
    	cout<<"Caught Exception: ParameterError";
    	return 0;
    }
    double r0=0,r1=0;
    //Already tested in EllipsePore
//    cout<<"Has cylindrical pore?\n";
//    cout<<"R up:";
//    cin>>r0;
//    cout<<"R bottom:";
//    cin>>r1;
    CylindricalHole hole=CylindricalHole(r0,r1,T+10*cutoff,-T/2-5*cutoff);
    if (r0 >0 && r1 >0) {
    pore.has_hole(&hole);
    }
	return 1;
}

int testConstructCharge(GridSpace& space,ToroidalPore& pore) {
	ChargeDensity c=ChargeDensity(2);
	c.set_nlayer(2);
	c.set_layer(0,3,0.75,-2.0*65);
	c.set_layer(1,1.5,0.75,1.0*65);
	double R0,K0,T;
	R0=pore.r0;
	K0=pore.k0;
	T=pore.T;
	double a,b;
	a=(20*R0*K0/T+6*K0*K0/T-5*R0-3*K0)/(10*R0+2*K0);
	b=4*(1-a)/(T*T);
	QuadraticChargeDistribution d;
	d.setup(a,b,R0,K0,T);
	pore.set_distribution(&c,&d);//,&d);
	GridSpace ref;
	clock_t t1=clock();
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Gaussian in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("add_charge.dx",&space);

	return (1);
}
int testConstructDiel(GridSpace& space,ToroidalPore& pore) {
	DielectricConstant c;
	//TODO: inside =6.0?
	c.set_diels(80.0,10.0,2.0,2.0);
	c.set_headgroup_width(3.0);
	double R0,K0,T;
	R0=pore.r0;
	K0=pore.k0;
	T=pore.T;
//	QuadraticDistribution d;
//	d.setup(0.2,1.3,R0,K0,T);
	pore.set_distribution(&c);
	GridSpace ref;
	clock_t t1=clock();
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Dielectric in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("add_diel.dx",&space);

	return (1);
}
int testConstructKappa(GridSpace& space,ToroidalPore& pore) {
	Kappa c;
	c.set_headgroup_width(3);
	double R0,K0,T;
	R0=pore.r0;
	K0=pore.k0;
	T=pore.T;
	QuadraticChargeDistribution d;
	d.setup(0.2,1.3,R0,K0,T);
	pore.set_distribution(&c);
	GridSpace ref;
	clock_t t1=clock();
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Kappa in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("add_kappa.dx",&space);
	return (1);
}

int testConstructIMM1(GridSpace& space,ToroidalPore& pore) {
	IMM1Transition c;
	double R0,K0,T;
	R0=pore.r0;
	K0=pore.k0;
	T=pore.T;
	c.set_nsmth(10);
	c.set_memb_thickness(T);
	QuadraticChargeDistribution d;
	d.setup(0.2,1.3,R0,K0,T);
	pore.set_distribution(&c);
	GridSpace ref;
	clock_t t1=clock();
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling IMM1 Transition function in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("IMM1.dx",&space);
	return (1);
}

int testConstructGeoConstraint1(GridSpace& space,ToroidalPore& pore) {
	GaussianGeoConstraint c;
	c.set_forceconstant(5.0);
	c.set_offset(3.0);
	c.set_width(5.0);
	double R0,K0,T;
	R0=pore.r0;
	K0=pore.k0;
	T=pore.T;
	pore.set_distribution(&c);
	GridSpace ref;
	clock_t t1=clock();
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Geo Constraint 1 in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("geo_constraint1.dx",&space);
	return (1);
}

int testConstructGeoConstraint2(GridSpace& space,ToroidalPore& pore) {
	SwitchGeoConstraint c;
	c.set_cutoff(10.0);
	c.set_forceconstant(5.0);
	double R0,K0,T;
	R0=pore.r0;
	K0=pore.k0;
	T=pore.T;
	QuadraticChargeDistribution d;
	d.setup(0.2,1.3,R0,K0,T);
	pore.set_distribution(&c);
	GridSpace ref;
	clock_t t1=clock();
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Geo Constraint 2 in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("geo_constraint2.dx",&space);
	return (1);
}

int testConstructGeoConstraint3(GridSpace& space,ToroidalPore& pore) {
	HarmonicGeoConstraint c;
	c.set_cutoff(3.0);
	c.set_forceconstant(5.0);
	c.set_offset(3.0);
	double R0,K0,T;
	R0=pore.r0;
	K0=pore.k0;
	T=pore.T;
	QuadraticChargeDistribution d;
	d.setup(0.2,1.3,R0,K0,T);
	pore.set_distribution(&c);
	GridSpace ref;
	clock_t t1=clock();
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Geo Constraint in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("geo_constraint3.dx",&space);
	return (1);
}

int main() {
	ToroidalPore pore;
	testSetupPore(pore);
	GridSpace space;
	DxFile f;
	DxFile::read_dx("charge_X.dx",&space);
	testConstructCharge(space,pore);
	f.load(&space);
	f.ave("charge.dat");
	DxFile::read_dx("dielx_X.dx",&space);
	testConstructDiel(space,pore);
	f.load(&space);
	f.ave("dielx.dat");
	DxFile::read_dx("kappa_X.dx",&space);
	testConstructKappa(space,pore);
	f.load(&space);
	f.ave("kappa.dat");
	space.clear(0);
	testConstructIMM1(space,pore);
	f.load(&space);
	f.ave("imm1.dat");
	space.clear(0);
	testConstructGeoConstraint1(space,pore);
	f.load(&space);
	f.ave("geo1.dat");
	space.clear(0);
	testConstructGeoConstraint2(space,pore);
	f.load(&space);
	f.ave("geo2.dat");
	space.clear(0);
	testConstructGeoConstraint3(space,pore);
	f.load(&space);
	f.ave("geo3.dat");
	return 0;
}
